International Cherry Blossom Prediction Competition

Exploratory Data Analysis 1

Document description

• En este documento se lleva a cabo análisis exploratorio de datos con las variables bioclimáticas, de suelo, coberturas terrestres y huella humana. Recordemos que las variables bioclimáticas son un consolidado de los años 1970 a 2000. Este documento tiene como finalidad comparar las fechas de floración entre países, validar relaciones que han sido reportadas en literatura (por ejemplo, la relación de temperatura con DOY), determinar la existencia de efectos de interacción y establecer cuáles de las variables predictoras podrían ser eventualmente candidatas para incluir en la modelación (selección de predictores). Podríamos decir que las variables analizadas en este documento son “estáticas” o “fijas”, ya que no tenemos temporalidad de ellas, son datos consolidados en ventanas de tiempo preestablecidas.
• En el siguiente documento (06-EDA2-FE.qmd) se muestran resultados de análisis exploratorio para variables “dinámicas”, me refiero a ellas como dinámicas ya que las variables climatológicas y de fotoperíodo tienen temporalidad (datos diarios) y obtuve información desde el 01 de enero de 1981 hasta el 31 de diciembre de 2023.

Libraries and setup

library(tidyverse)
library(arrow)
library(splines)
library(corrr)
library(scales)
library(FactoMineR)
library(factoextra)
library(glue)

theme_set(theme_bw() + theme(legend.position = "top"))

Data

• Extraigo una variable de nombre country que tiene el país, esto lo hago con la finalidad de evaluar patrones de comportamiento por diferentes países y porque podría ser de ayuda para extrapolar a regiones donde no tenemos registros históricos de floración. Para diferenciar los datos de NPN para USA, los etiqueto en la variable country como USA-NPN y a los de Washington DC como USA-WDC.
• En varios de los resultados de ese documento no se incluye USA-WDC porque sólo es una coordenada (lat=38.88535, long=-77.03863) y tampoco se muestran los resultados de USA-NPN, la razón es que las variables bioclimáticas no están para las coordenadas de USA-NPN, sin embargo, como información de clima fue suministrada por NPN, uso esta información para análisis exploratorio al final de este documento y también en el próximo (06-EDA2-FE.qmd).

df_full <- read_parquet("../external-data/df_full.parquet") |> 
  mutate(
    country = 
      case_when(
        str_detect(location, "Japan") ~ "Japan",
        str_detect(location, "kyoto") ~ "Japan",
        str_detect(location, "liestal") ~ "Switzerland",
        str_detect(location, "Switzerland") ~ "Switzerland",
        str_detect(location, "South Korea") ~ "South Korea",
        str_detect(location, "vancouver") ~ "Canada",
        str_detect(location, "washingtondc") ~ "USA-WDC",
        str_detect(location, "NPN") ~ "USA-NPN"
      )
  ) |> 
  relocate(location, country, everything())

DOY distribution by country

df_full |> 
  filter(location != "NPN") |> 
  ggplot(aes(x = bloom_doy)) +
  facet_wrap(~country, ncol = 1, scales = "free_y") +
  geom_histogram(color = "black") +
  geom_rug() +
  labs(x = "Bloom DOY", y = "Count",
       title = "DOY distribution by country",
       subtitle = "Original scale")
df_full |> 
  filter(location != "NPN") |> 
  ggplot(aes(x = bloom_doy)) +
  facet_wrap(~country, ncol = 1, scales = "free_y") +
  geom_histogram(color = "black") +
  geom_rug() +
  scale_x_log10() +
  labs(x = "Bloom DOY", y = "Count",
       title = "DOY distribution by country",
       subtitle = "Logarithmic scale")

Scatter plot of all vs DOY

• Los siguientes gráficos permiten evidenciar relaciones (lineales y no lineales) entre algunas de las variables bioclimáticas y el día de la floración, también se observan asociaciones entre composición de suelo, por ejemplo, pH y DOY, no son dicientes o claras las relaciones entre las variables de cobertura terreste, huella humana, pastizales, arbutos, agua y la fecha de la floración.
• Con estos gráficos se valida o confirma que la temperatura (bio1, bio5, bio8, bio9, bio10 y bio11) y el día de la floración exhiben una relación que podríamos sugerir está entre lineal y cuadrática, independiente del país, a mayor temperatura se han observado floraciones más precoces, este comportamiento es similar en Japón y Suiza, no obstante, Corea del Sur presenta algunas diferencias, por ejemplo, en la variable bio11 para Japón y Suiza, la relación de tipo cuadrático es evidente, pero en Corea del sur no lo es.
• Las precipitaciones (bio12 a bio19) también muestran relaciones interesantes con la fecha de floración, con ligeras diferencias entre países, por ejemplo, la precipitación anual (bio12) muestra una relación lineal inversa con el DOY, en Suiza por el contrario observamos una relación cuadrática positiva y en Corea del Sur no es clara la relación. Este resultado puede sugerir que temperatura y precipitación son factores que interactuan y tener un efecto conjunto; más adelante se exploran interacciones entre variables para determinar covariaciones de interés.
• La estacionalidad de la temperatura (bio4) muestra un patrón de comportamiento diferente entre países, en Japón la relación es positiva y lineal, en Corea del Sur no es clara la asociación y en Suiza muesra una relación inversa cuadrática. Recordemos que las variables bioclimáticas (bio_) recogen información de los años 1970 a 2000, con una resolución de \(1 km^2\), este hallazgo será comparado más adelante con la información diaria de temperatura para validar las relaciones que se evidencian en este gráfico.
• En cuanto a la composición del suelo, el nitrógeno, el pH del agua del suelo, el contenido de carbono orgánico del suelo (SOC), la densidad aparente (BDOD) y la capacidad de intercambio catiónico (CEC) muestran tendencias en la relación con el DOY. Las otras variables parecen no tener ninguna asociación (lineal o no lineal) con la variable objetivo.
• Finalmente, la información de tierras de cultivo (CROPLAND), huella humana (FOOTPRINT) y coberturas terrestres, no exhiben ningún tipo de relación que pudiera ser considerada de interés para análisis posteriores, además, algunas de estas variables (por ejemplo SHRUBS) carecen de información para las coordenadas bajo análisis.

Japan

df_full |>
  filter(country == "Japan") |> 
  select(-c(location, country,  year, bloom_date)) |>
  pivot_longer(cols = -bloom_doy) |>
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point(size = 0.5, alpha = 0.25) +
  stat_bin_2d(aes(fill = ..density..),
                  geom = "raster", 
                  contour = FALSE,
                  show.legend = FALSE) +
  scale_fill_distiller(palette = 4, direction = -1) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "firebrick3",
              se = FALSE) +
  scale_x_log10() +
  labs(y = "DOY", x = "", title = "Japan")

Switzerland

df_full |>
  filter(country == "Switzerland") |> 
  select(-c(location, country,  year, bloom_date)) |>
  pivot_longer(cols = -bloom_doy) |>
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point(size = 0.5, alpha = 0.25) +
  stat_bin_2d(aes(fill = ..density..),
                  geom = "raster", 
                  contour = FALSE,
                  show.legend = FALSE) +
  scale_fill_distiller(palette = 4, direction = -1) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "firebrick3",
              se = FALSE) +
  scale_x_log10() +
  labs(y = "DOY", x = "", title = "Switzerland")

South Korea

df_full |>
  filter(country == "South Korea") |>
  select(-c(location, country,  year, bloom_date)) |>
  pivot_longer(cols = -bloom_doy) |>
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point(size = 0.5, alpha = 0.25) +
  stat_bin_2d(aes(fill = ..density..),
                  geom = "raster", 
                  contour = FALSE,
                  show.legend = FALSE) +
  scale_fill_distiller(palette = 4, direction = -1) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "firebrick3",
              se = FALSE) +
  scale_x_log10() +
  labs(y = "DOY", x = "", title = "South Korea")

All the countries

df_full |>
  filter(location != "NPN") |>
  select(-c(location, country,  year, bloom_date)) |>
  pivot_longer(cols = -bloom_doy) |>
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point(size = 0.5, alpha = 0.25) +
  stat_bin_2d(aes(fill = ..density..),
                  geom = "raster", 
                  contour = FALSE,
                  show.legend = FALSE) +
  scale_fill_distiller(palette = 4, direction = -1) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "firebrick3",
              se = FALSE) +
  scale_x_log10() +
  labs(y = "DOY", x = "", title = "All the countries")

Non-parametric correlations

• Estos perfiles de correlación muestran para cada país y todos los países, de mayor a menor correlación (Spearman). En Japón y Corea del Sur, los cerezos ubicados más al norte parecen tardarse más en la floración (colores azules), la variable latitud es la primera en el perfil de correlaciones de estos países, sin embargo, en Suiza no es la latitud el factor lineal de mayor influcencia, es la altitud seguido de variables de composición de suelo y en el top 10 están tres variables que recogen información de precipitación. Si observamos la parte baja del gráfico (colores rojos) en todos los países coincide la variable temperatura anual (bio1) en ser el factor lineal inverso de mayor magnitud, seguido de otras variables indicadoras también de temperaturas. Este gráfico es de mucha utilidad no sólo para representar relaciones de tipo lineal sino también para obtener una caracterización rápida de la relación del clima con la fenología en cerezos ubicados en diferentes países, podemos decir que si pensamos en los factores que se asociación con tardanza en la floración, los perfiles entre países son muy diferentes, no obstante, si pensamos en los factores que se asociación con la precocidad en la floración, los perfiles entre países son muy parecidos. En conclusión, florecer más rápido parece tener unas “causas” comunes y globales ligadas a la temperatura, sin embargo, el retraso en la floración podría ser más difícil de establecer y generalizar.

df_full |>
  filter(country == "Japan") |>
  select(where(is.numeric), -c(shrubs)) |>
  correlate(method = "spearman") |>
  filter(term == "bloom_doy") |>
  pivot_longer(cols = -term) |>
  filter(!is.na(value)) |>
  mutate(
    name = str_replace_all(
      name,
      "wildareas-v3-1993-human-footprint",
      "human-footprint93"
    ),
    name = str_replace_all(
      name,
      "wildareas-v3-2009-human-footprint",
      "human-footprint09"
    )
  ) |>
  ggplot(aes(
    x = reorder(name, value),
    y = term,
    fill = value
  )) +
  geom_tile() +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1)) +
  labs(x = "",
       y = "",
       fill = "",
       title = "Japan") +
  coord_flip()
df_full |> 
  filter(country == "Switzerland") |> 
  select(where(is.numeric), -c(shrubs)) |> 
  correlate(method = "spearman") |> 
  filter(term == "bloom_doy") |> 
  pivot_longer(cols = -term) |> 
  filter(!is.na(value)) |> 
  mutate(name = str_replace_all(name,
                                "wildareas-v3-1993-human-footprint",
                                "human-footprint93"),
         name = str_replace_all(name,
                                "wildareas-v3-2009-human-footprint",
                                "human-footprint09")) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile() +
  geom_tile() +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1)) +
  labs(x = "", y = "", fill = "", title = "Switzerland") +
  coord_flip()
df_full |> 
  filter(country == "South Korea") |> 
  select(where(is.numeric), -c(shrubs)) |> 
  correlate(method = "spearman") |> 
  filter(term == "bloom_doy") |> 
  pivot_longer(cols = -term) |> 
  filter(!is.na(value)) |> 
  mutate(name = str_replace_all(name,
                                "wildareas-v3-1993-human-footprint",
                                "human-footprint93"),
         name = str_replace_all(name,
                                "wildareas-v3-2009-human-footprint",
                                "human-footprint09")) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile()  +
  geom_tile() +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1)) +
  labs(x = "", y = "", fill = "", title = "South Korea") +
  coord_flip()
df_full |> 
  filter(location != "NPN") |> 
  select(where(is.numeric), -c(shrubs)) |> 
  correlate(method = "spearman") |> 
  filter(term == "bloom_doy") |> 
  pivot_longer(cols = -term) |> 
  filter(!is.na(value)) |> 
  mutate(name = str_replace_all(name,
                                "wildareas-v3-1993-human-footprint",
                                "human-footprint93"),
         name = str_replace_all(name,
                                "wildareas-v3-2009-human-footprint",
                                "human-footprint09")) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile()  +
  geom_tile() +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1),
        legend.key.size = unit(0.85, "cm")) +
  labs(x = "", y = "", fill = "", title = "All the countries") +
  coord_flip()

Interaction of variables

• Con los siguientes gráficos pretendo reconocer si existen comportamientos de interacción, es decir, intento confirmar a partir de representaciones visuales si hay efectos conjuntos entre las variables predictoras. En vista de que el conjunto de variables predictoras es “grande” \((P = 41)\) generar todas las interacciones (dobles, triples o de mayor jerarquía) podría ser dispendioso y demorado. Las siguientes son las dos aproximaciones que implemento:
  • 1. Inicialmente voy a graficar diagramas de dispersión con las variables que mostraron tendencias en los gráficos de dispersión y los perfiles de correlación; el color de los puntos está determinado por la variable respuesta bloom_doy.
  • 2. Grafico la interacción doble para algunas variables de interés vs la variable respuesta.

Scatter plot

• En estos gráficos aplico la función distinct() por latitud y longitud porque la misma coordenada (así sean en años diferentes) tiene el mismo dato de bioclimáticas, suelo y coberturas, ya que esta información no posee temporalidad. Por esa razón resumo con la mediana del bloom_doy (ya que esta variable sí cambia con el tiempo).

Temperature - Precipitation

df_full |>
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Annual precipitation (mm)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Annual precipitation (mm)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Annual precipitation (mm)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Temperature - Nitrogen

df_full |>
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Temperature - SOC

df_full |>
  filter(`soc_0-5cm_mean` > 0) |> 
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Soil organic carbon (dg/kg)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`soc_0-5cm_mean` > 0) |> 
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Soil organic carbon (dg/kg)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`soc_0-5cm_mean` > 0) |> 
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Soil organic carbon (dg/kg)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Temperature - BDOD

df_full |>
  filter(`bdod_0-5cm_mean` > 0) |> 
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Bulk density (cg/cm³)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`bdod_0-5cm_mean` > 0) |> 
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Bulk density (cg/cm³)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`bdod_0-5cm_mean` > 0) |> 
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Bulk density (cg/cm³)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Precipitation - Nitrogen

df_full |>
  filter(`nitrogen_0-5cm_mean` > 0) |> 
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual precipitation (mm)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`nitrogen_0-5cm_mean` > 0) |> 
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual precipitation (mm)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`nitrogen_0-5cm_mean` > 0) |>  
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual precipitation (mm)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Altitude - Precipitation

df_full |>
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = alt, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(y = "Annual precipitation (mm)",
       x = "Altitude",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = alt, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(y = "Annual precipitation (mm)",
       x = "Altitude)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = alt, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(y = "Annual precipitation (mm)",
       x = "Altitude",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Interaction

En este caso como las variables son numéricas simplemento multiplico las dos variables bajo análisis para evaluar el efecto de interacción. Teniendo en cuenta que bloom_doy es la variable respuesta \((y)\) y si ejemplificamos un modelo sencillo que incluye las variables temperatura anual \((x_1)\) y precipitación \((x_2)\) como predictoras, mi intención es validar si el efecto acumulativo o conjunto está presente. Podemos estimar el efecto marginal de \(x_1\) y definirlo como \(\beta_1\), lo propio podemos hacer con el efecto marginal de \(x_2\) y definirlo como \(\beta_2\), sin embargo, también es de interés evaluar el efecto de interacción \((\beta_3)x_1x_2\), en ese orden de ideas se puede expresar el modelo matamático como se muestra a continuación:

\[y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2 + error\]

Un primer resultado interesante de las interaccione evaluadas es que el comportamiento discrepa entre países, por ejemplo, la interacción temperatura-precipitación muestra una relación inversa con bloom_doy, pero no es lineal en los tres países, en Corea del Sur la relación tiende a ser cuadrática. Algo similar podemos observar en la interacción temperatura-nitrógeno, en Japón y Corea del Sur no exhibe ninguna relación con la variable objetivo, sin embargo, en Suiza se evidencia una relación lineal inversa. En los tres países la interacción temperatura-bulk density se relaciona de forma negativa con la variable respuesta, pero en Japón es evidente la dependencia lineal que existe entre la variable \(y\) y la predictora derivada de esta interacción.

Más adelante (en el documento 06-EDA2-FE.qmd) hago uso de la regresión lasso para seleccionar predictores y evaluar posibles efectos de interacción que no fueron representados en este análisis exploratorio.

Temperature - Precipitation

countries <- c("Japan", "Switzerland", "South Korea")

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)",
         y = "DOY",
         title = countries[3])

Temperature - Nitrogen

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[3])

Temperature - SOC

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `soc_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `soc_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `soc_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)",
         y = "DOY",
         title = countries[3])

Temperature - BDOD

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)",
         y = "DOY",
         title = countries[3])

Precipitation - Nitrogen

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[3])

Altitude - Precipitation

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = alt * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Altitude * Annual precipitation (mm) ",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = alt * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Altitude * Annual precipitation (mm) ",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = alt * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Altitude * Annual precipitation (mm) ",
         y = "DOY",
         title = countries[3])

Principal Component Analysis (PCA)

• En estos gráficos aplico la función distinct() por latitud y longitud porque la misma coordenada (así sean en años diferentes) tiene el mismo dato de bioclimáticas, suelo y coberturas, ya que esta información no posee temporalidad. Por esa razón resumo con la mediana del bloom_doy (ya que esta variable sí cambia con el tiempo).
• Aprovechando que existe alta correlación entre las variables predictoras, ejecuto análisis de componentes principales para reducir la dimensionalidad e intentar ver si en un espacio más pequeño de variables se exhibe algún patrón de comportamiento relacionado con la variable respuesta.
• Para facilitar las visualizaciones discretizo la variable respuesta en los siguientes niveles (ordinales):
  • Q1: valores de bloom_doy inferiores al valor del cuartil 1: \(DOY < Q1\).
  • Q2: valores de bloom_doy mayores o iguales al valor del cuartil 1 y menores al valor del cuartil 2: \(Q1 \leq DOY < Q2\).
  • Q3: valores de bloom_doy mayores o iguales al valor del cuartil 2 y menores al valor del cuartil 3: \(Q2 \leq DOY < Q3\).
  • Q4: valores de bloom_doy superiores o iguales al valor del cuartil 3: \(DOY \geq Q3\).
• La variable respuesta discretizada doy_categ la introduzco al PCA como variable suplementaría cualitativa.
• Independientemente del país las floraciones que no superan el Q1 (DOY < Q1) y las que se tardan más del Q3 (DOY >= Q3) exhiben diferencias al graficarlas en los tres primeros componentes principales, quedando ubicadas en lugares diferentes. Es importante mencionar que el patrón es menos visible en Corea del Sur. Floraciones que están entre el el Q1 y Q3 tienden a ser similares. Esta diferenciación que se aprecia entre floraciones precoces y floraciones tardías podría constituirse en un indicativo de perfiles diferenciales en el nicho ecológico o entorno ambiental al cual están expuestos los cerezos.
• En el caso de Japón se observa que el componente 1 está asociado de forma positiva con variables que recogen información de temperatura (bio1, bio5, bio10, bio11) y precipitación (bio12, bio13, bio16), indicando que valores mayores a cero en el componente 1 los podemos asociar con valores altos de temperatura y precipitación, justamente bajo estas condiciones observamos las floraciones rápidas (puntos azules), de otro lado, este componente está asociado de forma negativa con la latitud, la longitud, bio4 (estacionalidad de la temperatura), bio7 (rango anual de temperatura) y en menor medida con el nitrógeno, indicando que valores mayores a cero de este componente estarán relacionados no sólo con valores altos de temperatura y precipitación sino también con menor variación en la temperatura (bio4) y menor rango de temperatura anual (bio7), por tanto, el perfil descrito es el que aplica a las floraciones precoces, justo lo contrario ocurre con floraciones tardías, donde se esperan temperaturas y precipitaciones bajas, con mayores variaciones en la temperatura a través del año, también podríamos afirmar que plantas que están ubicadas más al sur de Japón tienden a presentar floraciones retardadas. Por la ubicación de los puntos también podríamos afirmar que las floraciones precoces en Japón están moduladas por cambios de temperatura, sin embargo, otros factores como la capacidad de intercambio catiónico (cec) y la densidad del carbono orgánico (ocd) podrían ser influyentes, ya que se observa que estas variables tiene correlación positiva con el componente 2, indicando que valores mayores a cero en el componente 2 están asociados con suelos altos en capacidad de intercambio catiónico y suelos con mayor densidad del carbono orgánico.
• En Suiza también se evidencia la diferencia entre floraciones precoces y floraciones tardías, no obstante, la correlación de los componentes con las variables permite inferir patrones de comportamiento ligeramente diferentes a Japón y Corea del Sur. El componente uno tiene correlación positiva con variables que recogen información de temperatura (bio4, bio5, bio10, bio11) y exhibe correlación negativa con la altitud, el contenido de carbono orgánico del suelo (soc), el inventario de carbono orgánico del suelo (ocs) y en menor medida con variables de precipitación (bio12, bio16, bio17), estas correlaciones permiten inferior que floraciones precoses se dan en condiciones de alta temperatura (también la variación anual) y en sitios donde hay menor cantidad de lluvia anual, menor contenido de carbono orgánico del suelo y ubicadas en altitudes bajas. Los cerezos que exhiben floraciones tardías tienden a estar ubicados en altitudes elevadas, sitios con mayor cantidad de precipitación anual, mayores contenidos de carbono orgánico del suelo.
• En Corea del Sur se observa un comportamiento similar al de Japón y Suiza en cuanto a que la temperatura es el factor determinante en los tiempos de floración. Sin embargo, el componente 1 está asociado de forma negativa con variables que recogen información de temperatura (bio1, bio6, bio9, bio11) y de forma positiva con variables que recogen información de la variación anual de la temperatura (bio4, bio7). Las floraciones tardías se ubican a la derecha del 0 en el componente 1, es decir, que estas plantas se espera estén en lugares con mayores variaciones en la temperatura a través del año y también se espera que las temperaturas de estos sitios sean más bajas. Basado en la asociación positiva que tiene el componente 3 con las precipitaciones (bio12, bio15, bio16) y observando la ubicación de los puntos verdes en el gráfico de CP1 vs CP3, es posible intuir que las floraciones demoradas también están asociadas con mayores precipitaciones.
• Es importante mencionar que todas estas interpretaciones las hago desde la perspectiva exploratoria, siguiendo las premisas de J. W. Tukey como se muestra en el siguiente diagrama. La idea principal ha sido reconocer patrones que me permitan realizar un diagnóstico objetivo de cara a la fase de modelación.

Japan

• Con cuatro componentes principales se retiene el 69.4% de la variabilidad total.
• Con diez componentes principales se retiene el 89.2% de la variabilidad total.
• Disminuir de 41 dimensiones a 4 perdiente aproximadamente el 30% de la variabilidad no está nada mal para propósitos de graficación, sin embargo, para propósitos de modelación, valdría la pena reducir de 41 variables a solo 10 componentes.

df_pca_japan1 <-
  df_full |>
  filter(country == countries[1]) 

value_q1 <- 
  quantile(df_pca_japan1$bloom_doy, probs = 0.25, na.rm = TRUE)

value_q2 <- 
  quantile(df_pca_japan1$bloom_doy, probs = 0.50, na.rm = TRUE)

value_q3 <- 
  quantile(df_pca_japan1$bloom_doy, probs = 0.75, na.rm = TRUE)

order_categ <-
  c("DOY < Q1",
    "Q1 <= DOY < Q2",
    "Q2 <= DOY < Q3",
    "DOY >= Q3")

df_pca_japan2 <-
  df_pca_japan1 |>
  group_by(lat, long) |>
  mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
  ungroup() |>
  distinct(lat, long, .keep_all = TRUE) |>
  select(-shrubs) |>
  mutate(
    doy_categ = case_when(
      median_doy < value_q1 ~ "DOY < Q1",
      median_doy >= value_q1 &
        median_doy < value_q2 ~ "Q1 <= DOY < Q2",
      median_doy >= value_q2 &
        median_doy < value_q3 ~ "Q2 <= DOY < Q3",
      median_doy >= value_q3 ~ "DOY >= Q3",
    ),
    doy_categ = factor(doy_categ, levels = order_categ)
  ) |>
  select(-c(location, country, bloom_date, bloom_doy, median_doy)) |>
  relocate(doy_categ, everything())


pca_japan <- PCA(X = df_pca_japan2, graph = FALSE, quali.sup = 1)

df_pca_japan2$pc1 <- pca_japan$ind$coord[, 1]
df_pca_japan2$pc2 <- pca_japan$ind$coord[, 2]
df_pca_japan2$pc3 <- pca_japan$ind$coord[, 3]

eigen_pc1 <- pca_japan$eig[, 2][1] |> round(digits = 1)
eigen_pc2 <- pca_japan$eig[, 2][2] |> round(digits = 1)
eigen_pc3 <- pca_japan$eig[, 2][3] |> round(digits = 1)

df_pca_japan2 |> 
  ggplot(aes(x = pc1, y = pc2, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 2 ({eigen_pc2}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 2")
df_pca_japan2 |>
  select(-c(doy_categ)) |> 
  correlate(method = "pearson") |> 
  filter(term %in% c("pc1", "pc2", "pc3")) |> 
  select(-c(pc1, pc2, pc3)) |> 
  pivot_longer(-term) |> 
  mutate(
    name = str_replace_all(
      name,
      "wildareas-v3-1993-human-footprint",
      "human-footprint93"
    ),
    name = str_replace_all(
      name,
      "wildareas-v3-2009-human-footprint",
      "human-footprint09"
    )
  ) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile() + 
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  labs(x = "", y = "Component", fill = "Correlation") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5),
        legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))
df_pca_japan2 |> 
  ggplot(aes(x = pc1, y = pc3, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 3 ({eigen_pc3}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 3")

Switzerland

• Con cuatro componentes principales se retiene el 66.4% de la variabilidad total.
• Con diez componentes principales se retiene el 85.5% de la variabilidad total.

df_pca_switzerland1 <-
  df_full |>
  filter(country == countries[2]) 

value_q1 <- 
  quantile(df_pca_switzerland1$bloom_doy, probs = 0.25, na.rm = TRUE)

value_q2 <- 
  quantile(df_pca_switzerland1$bloom_doy, probs = 0.50, na.rm = TRUE)

value_q3 <- 
  quantile(df_pca_switzerland1$bloom_doy, probs = 0.75, na.rm = TRUE)

order_categ <-
  c("DOY < Q1",
    "Q1 <= DOY < Q2",
    "Q2 <= DOY < Q3",
    "DOY >= Q3")

df_pca_switzerland2 <-
  df_pca_switzerland1 |>
  group_by(lat, long) |>
  mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
  ungroup() |>
  distinct(lat, long, .keep_all = TRUE) |>
  select(-shrubs) |>
  mutate(
    doy_categ = case_when(
      median_doy < value_q1 ~ "DOY < Q1",
      median_doy >= value_q1 &
        median_doy < value_q2 ~ "Q1 <= DOY < Q2",
      median_doy >= value_q2 &
        median_doy < value_q3 ~ "Q2 <= DOY < Q3",
      median_doy >= value_q3 ~ "DOY >= Q3",
    ),
    doy_categ = factor(doy_categ, levels = order_categ)
  ) |>
  select(-c(location, country, bloom_date, bloom_doy, median_doy)) |>
  relocate(doy_categ, everything())


pca_switzerland <- PCA(X = df_pca_switzerland2, graph = FALSE, quali.sup = 1)

df_pca_switzerland2$pc1 <- pca_switzerland$ind$coord[, 1]
df_pca_switzerland2$pc2 <- pca_switzerland$ind$coord[, 2]
df_pca_switzerland2$pc3 <- pca_switzerland$ind$coord[, 3]

eigen_pc1 <- pca_switzerland$eig[, 2][1] |> round(digits = 1)
eigen_pc2 <- pca_switzerland$eig[, 2][2] |> round(digits = 1)
eigen_pc3 <- pca_switzerland$eig[, 2][3] |> round(digits = 1)

df_pca_switzerland2 |> 
  ggplot(aes(x = pc1, y = pc2, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 2 ({eigen_pc2}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 2")
df_pca_switzerland2 |>
  select(-c(doy_categ)) |> 
  correlate(method = "pearson") |> 
  filter(term %in% c("pc1", "pc2", "pc3")) |> 
  select(-c(pc1, pc2, pc3)) |> 
  pivot_longer(-term) |> 
  mutate(
    name = str_replace_all(
      name,
      "wildareas-v3-1993-human-footprint",
      "human-footprint93"
    ),
    name = str_replace_all(
      name,
      "wildareas-v3-2009-human-footprint",
      "human-footprint09"
    )
  ) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile() + 
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  labs(x = "", y = "Component", fill = "Correlation") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5),
        legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))
df_pca_switzerland2 |> 
  ggplot(aes(x = pc1, y = pc3, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 3 ({eigen_pc3}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 3")

South Korea

• Con cuatro componentes principales se retiene el 66.4% de la variabilidad total.
• Con diez componentes principales se retiene el 89.7% de la variabilidad total.

df_pca_southk1 <-
  df_full |>
  filter(country == countries[3]) 

value_q1 <- 
  quantile(df_pca_southk1$bloom_doy, probs = 0.25, na.rm = TRUE)

value_q2 <- 
  quantile(df_pca_southk1$bloom_doy, probs = 0.50, na.rm = TRUE)

value_q3 <- 
  quantile(df_pca_southk1$bloom_doy, probs = 0.75, na.rm = TRUE)

order_categ <-
  c("DOY < Q1",
    "Q1 <= DOY < Q2",
    "Q2 <= DOY < Q3",
    "DOY >= Q3")

df_pca_southk2 <-
  df_pca_southk1 |>
  group_by(lat, long) |>
  mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
  ungroup() |>
  distinct(lat, long, .keep_all = TRUE) |>
  select(-shrubs) |>
  mutate(
    doy_categ = case_when(
      median_doy < value_q1 ~ "DOY < Q1",
      median_doy >= value_q1 &
        median_doy < value_q2 ~ "Q1 <= DOY < Q2",
      median_doy >= value_q2 &
        median_doy < value_q3 ~ "Q2 <= DOY < Q3",
      median_doy >= value_q3 ~ "DOY >= Q3",
    ),
    doy_categ = factor(doy_categ, levels = order_categ)
  ) |>
  select(-c(location, country, bloom_date, bloom_doy, median_doy)) |>
  relocate(doy_categ, everything())


pca_southk <- PCA(X = df_pca_southk2, graph = FALSE, quali.sup = 1)

df_pca_southk2$pc1 <- pca_southk$ind$coord[, 1]
df_pca_southk2$pc2 <- pca_southk$ind$coord[, 2]
df_pca_southk2$pc3 <- pca_southk$ind$coord[, 3]

eigen_pc1 <- pca_southk$eig[, 2][1] |> round(digits = 1)
eigen_pc2 <- pca_southk$eig[, 2][2] |> round(digits = 1)
eigen_pc3 <- pca_southk$eig[, 2][3] |> round(digits = 1)

df_pca_southk2 |> 
  ggplot(aes(x = pc1, y = pc2, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 2 ({eigen_pc2}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 2")
df_pca_southk2 |>
  select(-c(doy_categ)) |> 
  correlate(method = "pearson") |> 
  filter(term %in% c("pc1", "pc2", "pc3")) |> 
  select(-c(pc1, pc2, pc3)) |> 
  pivot_longer(-term) |> 
  mutate(
    name = str_replace_all(
      name,
      "wildareas-v3-1993-human-footprint",
      "human-footprint93"
    ),
    name = str_replace_all(
      name,
      "wildareas-v3-2009-human-footprint",
      "human-footprint09"
    )
  ) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile() + 
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  labs(x = "", y = "Component", fill = "Correlation") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5),
        legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))
df_pca_southk2 |> 
  ggplot(aes(x = pc1, y = pc3, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 3 ({eigen_pc3}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 3")

USA National Phenology Network

• En este punto voy a realizar análisis exploratorio con los datos suministrados por USA National Phenology Network.
• Como no hay información de variables bioclimáticas para estas coordenadas voy a agregar la información climática que fue suministrada inicialmente. Tampoco se tiene en cuenta la variable shrubs porque todos sus valores están en cero.
• Desconozco el significado o implicación práctica de algunas variables de los datos iniciales suministrados por NPN, también hay variables que no son necesarias para el abordaje que implemento, por tal motivo no tengo en cuenta las siguientes variables:
  • Phenophase_Description
  • Species_ID
  • Genus
  • Common_Name
  • Kingdom
  • Phenophase_ID
  • Phenophase_Description
  • NumDays_Since_Prior_No
  • First_Yes_Year
  • First_Yes_Month
  • First_Yes_Day
  • First_Yes_Julian_Date
  • Last_Yes_Julian_Date
  • Last_Yes_Year
  • Last_Yes_Month
  • Last_Yes_Day
  • Last_Yes_Julian_Date
  • NumDays_Until_Next_No
    
  • AGDD_in_F
  • Observation_Date
  • Last_Yes_DOY
  • First_Yes_DOY
• Con la finalidad de unificar la variable objetivo bloom_doy utilizo el promedio de las variables First_Yes_DOY y Last_Yes_DOY en los datos USA-NPN_individual_phenometrics_data.csv. Estas dos variables también hacen parte de la lista de variables excluidas.
• Uno las dos bases de datos suministradas por NPN y a esta información le agrego los datos fijos que se obtuvieron para cada coordenada. La base de datos final queda con 4101 filas (en total 152 coordenadas diferentes) y 43 columnas.
• Con ayuda de Google Maps y bboxfinder Verifiqué cada una de las 152 coordenadas únicas y extraje dos variables a partir de esas coordenadas: el país y la ciudad (contiene el estado). Esta información la agrego a los datos iniciales para realizar comparaciones por países o estados de un mismo país.

Data integration

# External data NPN
data_npn_fixed <-
  df_full |>
  filter(country == "USA-NPN") |>
  select(-contains("bio"),
         -c(shrubs, location, country, year, bloom_date,
            bloom_doy, alt)) |> 
  distinct(lat, long, .keep_all = TRUE)

# NPN
var_drop <-
  c(
    "Phenophase_Description",
    "Species_ID",
    "Genus",
    "Common_Name",
    "Kingdom",
    "Phenophase_ID",
    "Phenophase_Description",
    "NumDays_Since_Prior_No",
    "First_Yes_Year",
    "First_Yes_Month",
    "First_Yes_Day",
    "First_Yes_Julian_Date",
    "Last_Yes_Julian_Date",
    "Last_Yes_Year",
    "Last_Yes_Month",
    "Last_Yes_Day",
    "Last_Yes_Julian_Date",
    "NumDays_Until_Next_No",
    "AGDD_in_F",
    "Last_Yes_DOY",
    "First_Yes_DOY"
  )

data_npn1 <-
  read_csv("../data/USA-NPN_individual_phenometrics_data.csv") |>
  mutate(
    Observation_Date = str_c(First_Yes_Month, "-", First_Yes_Day, "-", First_Yes_Year),
    Observation_Date = mdy(Observation_Date),
    bloom_doy = (Last_Yes_DOY + First_Yes_DOY) / 2
  ) |>
  distinct_all() |>
  rename(lat = Latitude,
         long = Longitude,
         alt = Elevation_in_Meters) |>
  select(!all_of(var_drop)) |> 
  relocate(Observation_Date,
           Site_ID,
           Individual_ID,
           Species,
           bloom_doy,
           everything())

data_npn2 <- 
  read_csv("../data/USA-NPN_status_intensity_observations_data.csv") |> 
  mutate(
    Observation_Date = mdy(Observation_Date),
    bloom_doy = Day_of_Year
  ) |>
  distinct_all() |>
  rename(lat = Latitude,
         long = Longitude,
         alt = Elevation_in_Meters) |>
  select(all_of(names(data_npn1))) |> 
  relocate(Observation_Date,
           Site_ID,
           Individual_ID,
           Species,
           bloom_doy,
           everything())

# Join data NPN + data "fixed" or "static"
data_npn <- bind_rows(data_npn1, data_npn2) |> 
  left_join(data_npn_fixed, by = c("lat", "long")) |> 
  distinct_all()

# Countries vector
countries <-
  c(
    "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "Canada",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "Canada",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "China",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "Japan",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "Canada",
     "USA",
     "USA",
     "USA",
     "USA",
     "Canada",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "China",
     "USA",
     "USA",
     "China",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "USA",
     "Canada",
     "USA",
     "USA",
     "USA",
     "USA",
     "Canada",
     "Canada",
     "USA"
  )

# Cities vector
cities <-
  c(
    "Fayetteville, North Carolina",
     "Raleigh, North Carolina",
     "Washington D.C.",
     "Knoxville, Tennessee",
     "San francisco California",
     "Roanoke, Virginia",
     "Burlington, Vermont",
     "Gatlinburg, Tennessee",
     "Pigeon Forge, Tennessee",
     "Denver, Colorado",
     "New York, New York",
     "Evansville, Indiana",
     "Boston, Massachusetts",
     "Fargo, North Dakota",
     "Charlotte, North Carolina",
     "Kansas City, Missouri",
     "New York, New York",
     "Duluth, Minnesota",
     "Thunder Bay, Ontario",
     "Washington D.C.",
     "Washington D.C.",
     "Washington D.C.",
     "Washington D.C.",
     "Washington D.C.",
     "Newport, Oregon",
     "Washington D.C.",
     "Washington D.C.",
     "Washington D.C.",
     "Washington D.C.",
     "Washington D.C.",
     "Washington D.C.",
     "Washington D.C.",
     "Alexandria, Virginia",
     "Alexandria, Virginia",
     "New York, New York",
     "New York, New York",
     "New York, New York",
     "Alexandria, Virginia",
     "New York, New York",
     "Minneapolis, Minnesota",
     "London, Kentucky",
     "Hagerstown, Maryland",
     "London, Kentucky",
     "Nashville, Tennessee",
     "Philadelphia, Pennsylvania",
     "Allentown, Pennsylvania",
     "New York, New York",
     "Sudbury, Ontario",
     "Philadelphia, Pennsylvania",
     "Roanoke, Virginia",
     "Denver, Colorado",
     "Lexington, Kentucky",
     "Indianapolis, Indiana",
     "Denver, Colorado",
     "Chengdu, Sichuan",
     "New York, New York",
     "Sacramento, California",
     "Asheville, North Carolina",
     "New York, New York",
     "Salem, Oregon",
     "Seattle, Washington",
     "Seattle, Washington",
     "Salem, Oregon",
     "Portland, Oregon",
     "Seattle, Washington",
     "Seattle, Washington",
     "Tacoma, Washington",
     "Muncie, Ind.",
     "Kyoto, Kyoto",
     "Redmond, Washington",
     "Salem, Oregon",
     "Northampton, Massachusetts",
     "Asheville, North Carolina",
     "Salem, Oregon",
     "Raleigh, North Carolina",
     "Albany, New York",
     "Duluth, Minnesota",
     "New York, New York",
     "New York, New York",
     "Providence, Rhode Island",
     "Seattle, Washington",
     "Greeley, Colorado",
     "Lexington, Kentucky",
     "Portland, Oregon",
     "Columbus, Ohio",
     "North Bend, Washington",
     "Denver, Colorado",
     "Seattle, Washington",
     "Sacramento, California",
     "Klamath Falls, Oregon",
     "Reno, Nevada",
     "Edmonton, Alberta",
     "Eugene, Oregon",
     "Corvallis, Oregon",
     "Longview, Washington",
     "Port Angeles, Washington",
     "Vancouver, British Columbia",
     "Salem, Oregon",
     "Albany, Oregon",
     "Portland, Oregon",
     "Lebanon, Virginia",
     "Asheville, North Carolina",
     "Boston, Massachusetts",
     "Baltimore, Maryland",
     "Newark, New Jersey",
     "San francisco California",
     "New York, New York",
     "Okinawa",
     "Lexington, Kentucky",
     "New York, New York",
     "Dalian, Liaoning",
     "Fayetteville, North Carolina",
     "Raleigh, North Carolina",
     "Goldsboro, North Carolina",
     "Buffalo, New York",
     "Newburgh, New York",
     "Amherst, Massachusetts",
     "Buffalo, New York",
     "New York, New York",
     "Washington D.C.",
     "Providence, Rhode Island",
     "Cambridge, Massachusetts",
     "Atlantic City, New Jersey",
     "Newport, Oregon",
     "Lexington, Kentucky",
     "New York, New York",
     "Portland, Oregon",
     "Greenfield, Massachusetts",
     "Riverhead, New York",
     "Raleigh, North Carolina",
     "New York, New York",
     "New York, New York",
     "Allentown, Pennsylvania",
     "Williamsport, Pennsylvania",
     "Baltimore, Maryland",
     "Salt Lake City, Utah",
     "Lexington, Kentucky",
     "Providence, Rhode Island",
     "Newark, New Jersey",
     "Harrisburg, Pennsylvania",
     "Philadelphia, Pennsylvania",
     "Johnson City, Tennessee",
     "Syracuse, New York",
     "Wilmington, Delaware",
     "Vancouver, British Columbia",
     "Chicago, Illinois",
     "Burlington, Vermont",
     "Boston, Massachusetts",
     "Portland, Oregon",
     "Toronto, Ontario",
     "Niagara Falls, Ontario",
     "Eugene, Oregon"
  )

data_npn <-
  data_npn |>
  left_join(
    data_npn |>
      distinct(State, lat, long) |>
      mutate(country = countries,
             city = cities),
    by = c("State", "lat", "long")
  )

# Export data
write_parquet(data_npn, "../external-data/df_npn_usa.parquet")

data_npn |> head()

Distribution of variables (external data)

• Nota: como estas variables son fijas, para cada coordenada se filtraron registros diferentes para latitud y longitud.

data_npn_fixed |> 
  pivot_longer(cols = where(is.numeric)) |> 
  ggplot(aes(x = value)) +
  facet_wrap(~name, scales = "free", ncol = 4) +
  geom_histogram(color = "black")

Distribution of variables (internal data)

• Para este gráfico filtro los datos que sean mayores a -100, ya que por alguna razón hay datos con valores negativos. Supongo que son registros con ausencia de datos y por esa razón se asingó dicho valor (-9999).
• Estas variables son dinámicas y cambian con el tiempo, a diferencia de las anteriores que son fijas.

data_npn |> 
  select(!all_of(names(data_npn_fixed)))|> 
  select(-c(Site_ID, Individual_ID)) |> 
  pivot_longer(cols = where(is.numeric)) |> 
  filter(value > -100) |> 
  ggplot(aes(x = value)) +
  facet_wrap(~name, scales = "free", ncol = 4) +
  geom_histogram(color = "black")

DOY distribution by species

• En el gráfico anterior se observa alta dispersión en la variable respuesta bloom_doy, por esa razón en el siguiente gráfico comparo la distribución por especie.

data_npn |> 
  ggplot(aes(x = bloom_doy)) +
  facet_wrap(~Species, ncol = 1, scales = "free_y") +
  geom_histogram(color = "black") +
  geom_rug() +
  labs(x = "Bloom DOY", y = "Count",
       title = "DOY distribution by specie",
       subtitle = "Original scale")
data_npn |> 
  ggplot(aes(x = bloom_doy)) +
  facet_wrap(~Species, ncol = 1, scales = "free_y") +
  geom_histogram(color = "black") +
  geom_rug() +
  scale_y_log10() +
  labs(x = "Bloom DOY", y = "Count",
       title = "DOY distribution by specie",
       subtitle = "Original scale")

DOY distribution by species and state

• En el gráfico anterior se sigue observando alta dispersión en la variable respuesta bloom_doy, por esa razón en el siguiente gráfico comparo la distribución por especie y State.
• Filtro estados con 5 o más registros para contar con un mínimo de información y lograr representar la distribución.

list_species <- unique(data_npn$Species)

yedoensis

data_npn |>
  filter(Species == list_species[1]) |>
  mutate(State = str_c(State, " - ", city)) |>
  group_by(State) |>
  mutate(flag_n = n()) |>
  ungroup() |>
  filter(flag_n >= 5) |>
  ggplot(aes(x = State, y = bloom_doy)) +
  geom_boxplot() +
  labs(
    y = "Bloom DOY",
    x = "",
    title = "DOY distribution by state",
    subtitle = glue("{list_species[1]}")
  ) +
  coord_flip() +
  scale_y_log10()

serrulata

data_npn |>
  filter(Species == list_species[2]) |>
  mutate(State = str_c(State, " - ", city)) |>
  group_by(State) |>
  mutate(flag_n = n()) |>
  ungroup() |>
  filter(flag_n >= 5) |>
  ggplot(aes(x = State, y = bloom_doy)) +
  geom_boxplot() +
  labs(
    y = "Bloom DOY",
    x = "",
    title = "DOY distribution by state",
    subtitle = glue("{list_species[2]}")
  ) +
  coord_flip() +
  scale_y_log10()

pensylvanica

data_npn |>
  filter(Species == list_species[3]) |>
  mutate(State = str_c(State, " - ", city)) |>
  group_by(State) |>
  mutate(flag_n = n()) |>
  ungroup() |>
  filter(flag_n >= 5) |>
  ggplot(aes(x = State, y = bloom_doy)) +
  geom_boxplot() +
  labs(
    y = "Bloom DOY",
    x = "",
    title = "DOY distribution by state",
    subtitle = glue("{list_species[3]}")
  ) +
  coord_flip() +
  scale_y_log10()

cerasus

data_npn |>
  filter(Species == list_species[4]) |>
  mutate(State = str_c(State, " - ", city)) |>
  group_by(State) |>
  mutate(flag_n = n()) |>
  ungroup() |>
  filter(flag_n >= 5) |>
  ggplot(aes(x = State, y = bloom_doy)) +
  geom_boxplot() +
  labs(
    y = "Bloom DOY",
    x = "",
    title = "DOY distribution by state",
    subtitle = glue("{list_species[4]}")
  ) +
  coord_flip() +
  scale_y_log10()

emarginata

data_npn |>
  filter(Species == list_species[5]) |>
  mutate(State = str_c(State, " - ", city)) |>
  group_by(State) |>
  mutate(flag_n = n()) |>
  ungroup() |>
  filter(flag_n >= 5) |>
  ggplot(aes(x = State, y = bloom_doy)) +
  geom_boxplot() +
  labs(
    y = "Bloom DOY",
    x = "",
    title = "DOY distribution by state",
    subtitle = glue("{list_species[5]}")
  ) +
  coord_flip() +
  scale_y_log10()

Temperature vs DOY

yedoensis

data_npn |> 
  filter(Species == list_species[1]) |>
  select(bloom_doy, contains("Tm")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[1]}")) +
  scale_y_log10() 

serrulata

data_npn |> 
  filter(Species == list_species[2]) |>
  select(bloom_doy, contains("Tm")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[2]}")) +
  scale_y_log10() 

pensylvanica

data_npn |> 
  filter(Species == list_species[3]) |>
  select(bloom_doy, contains("Tm")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[3]}")) +
  scale_y_log10() 

cerasus

data_npn |> 
  filter(Species == list_species[4]) |>
  select(bloom_doy, contains("Tm")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[4]}")) +
  scale_y_log10() 

emarginata

data_npn |> 
  filter(Species == list_species[5]) |>
  select(bloom_doy, contains("Tm")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[5]}")) +
  scale_y_log10() 

Precipitation vs DOY

yedoensis

data_npn |> 
  filter(Species == list_species[1]) |>
  select(bloom_doy, contains("Prcp_")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[1]}")) +
  scale_y_log10() 

serrulata

data_npn |> 
  filter(Species == list_species[2]) |>
  select(bloom_doy, contains("Prcp_")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[2]}")) +
  scale_y_log10() 

pensylvanica

data_npn |> 
  filter(Species == list_species[3]) |>
  select(bloom_doy, contains("Prcp_")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[3]}")) +
  scale_y_log10() 

cerasus

data_npn |> 
  filter(Species == list_species[4]) |>
  select(bloom_doy, contains("Prcp_")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[4]}")) +
  scale_y_log10() 

emarginata

data_npn |> 
  filter(Species == list_species[5]) |>
  select(bloom_doy, contains("Prcp_")) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[5]}")) +
  scale_y_log10() 

Others vs DOY

• La asimetría en estas variables dificulta o

yedoensis

data_npn |> 
  filter(Species == list_species[1]) |>
  select(bloom_doy, lat, long, alt, AGDD, Daylength) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 5) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[1]}")) +
  scale_y_log10() 

serrulata

data_npn |> 
  filter(Species == list_species[2]) |>
  select(bloom_doy, lat, long, alt, AGDD, Daylength) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 5) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[2]}")) +
  scale_y_log10() 

pensylvanica

data_npn |> 
  filter(Species == list_species[3]) |>
  select(bloom_doy, lat, long, alt, AGDD, Daylength) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 5) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[3]}")) +
  scale_y_log10() 

cerasus

data_npn |> 
  filter(Species == list_species[4]) |>
  select(bloom_doy, lat, long, alt, AGDD, Daylength) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 5) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[4]}")) +
  scale_y_log10() 

emarginata

data_npn |> 
  filter(Species == list_species[5]) |>
  select(bloom_doy, lat, long, alt, AGDD, Daylength) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 5) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[5]}")) +
  scale_y_log10() 

Fixed variables vs DOY

yedoensis

data_npn |> 
  filter(Species == list_species[1]) |>
  select(bloom_doy, `bdod_0-5cm_mean`:water) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[1]}")) +
  scale_y_log10()

serrulata

data_npn |> 
  filter(Species == list_species[2]) |>
  select(bloom_doy, `bdod_0-5cm_mean`:water) |>  
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[2]}")) +
  scale_y_log10() 

pensylvanica

data_npn |> 
  filter(Species == list_species[3]) |>
  select(bloom_doy, `bdod_0-5cm_mean`:water) |> 
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[3]}")) +
  scale_y_log10() 

cerasus

data_npn |> 
  filter(Species == list_species[4]) |>
  select(bloom_doy, `bdod_0-5cm_mean`:water) |>  
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[4]}")) +
  scale_y_log10() 

emarginata

data_npn |> 
  filter(Species == list_species[5]) |>
  select(bloom_doy, `bdod_0-5cm_mean`:water) |>  
  pivot_longer(cols = -c(bloom_doy))  |>
  filter(value > -100) |> 
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point()  +
  geom_smooth(
    method = "gam",
    formula = y ~ ns(x, df = 2),
    color = "red",
    size = 0.5
  ) +
  labs(x = "", y = "DOY",
       title = glue("{list_species[5]}")) +
  scale_y_log10()